############################################################################################################################
########Inclusion, recognition, and inter-group comparisons: The effects of power-sharing institutions on grievances########
#####Individual-level analyses - predicted probability plots#####
#####Andreas Juon, 2023 (Journal of Conflict Research)#####
############################################################################################################################

###note: run this file after running "juon_jcr_grievances_analysis_ind.R"

####################################################################################
#########STEP 1: LIBRARIES, DATA IMPORT, DEFINITIONS#########
####################################################################################

setwd("XXXXX")#set to directory that contains the data files

#######1.1 Libraries#######
library(lme4)
library(boot)
library(plyr)
library(dplyr)
library(broom)

#######1.2 Defining variable lists#######
categorical_vars <- c("state_control","included_nc","ongoing_recent_contestation","educ_high","interest_pol_d")
nominal_vars <- c("region","survey","gender","cowcode","cowcode_year","gwgroupid_year","other")
numeric_vars <- c("age","ps_corp","ps_corp_all_diff_neg","tek_ps_corp_diff_neg","ps_lib","size","cpi","politya","lgdppc","fractionalization1")

######1.3 Formula for combining ggplot legends######
grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {
  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + 
                    theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x +
                 theme(legend.position = "none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)
  
  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl), 
                                            legend,ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend, ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))
  
  grid.newpage()
  grid.draw(combined)
  
  invisible(combined)
}

####################################################################################
#########STEP 2: Predicted probability plots (government dissatisfaction)#########
####################################################################################

######2.1 mean/median/modal values, all/included groups, used for predicted probabilities######
unsat_gov_means <- subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50)
unsat_gov_means <- subset(unsat_gov_means, state_control == 0)
unsat_gov_means <- unsat_gov_means[c("survey","cowcode","cowcode_year","gwgroupid_year","state_control","included_nc","ps_corp","ps_corp_all_diff_neg","tek_ps_corp_diff_neg","ps_lib",controls_gc,controls_i)]
for (i in (categorical_vars)) {
  step1 <- paste("unsat_gov_means$",i,"<- median(unsat_gov_means$",i,")",sep="")
  eval(parse(text=c(step1)))
}
for (i in (nominal_vars)) {
  step1 <- paste("unsat_gov_means$",i,"<- mlv(unsat_gov_means$",i,")",sep="")
  eval(parse(text=c(step1)))
}
for (i in (numeric_vars)) {
  step1 <- paste("unsat_gov_means$",i,"<- median(unsat_gov_means$",i,")",sep="")
  eval(parse(text=c(step1)))
}
unsat_gov_means <- unique(unsat_gov_means)
unsat_gov_means$merge <- 1

######2.2 dataframe: scenarios######
#for m1_unsat_gov: included_nc
m1_unsat_gov_pred_df <- left_join(data.frame(merge = 1, outcome = "dissatisfied with government", type ="domestic", included_nc = c(0, 0, 1), state_control = c(1,0,0), group_type = c("core","excluded NC","included NC")), subset(unsat_gov_means, select=-c(included_nc, state_control)), by="merge")
m1_unsat_gov_pred_df <- unique(m1_unsat_gov_pred_df)
#for m2_unsat_gov: ps_corp
m2_unsat_gov_pred_df <- left_join(data.frame(merge = 1, outcome = "dissatisfied with government", type ="domestic", ps_corp = rep(seq(0, 1, length.out = 50), 2), state_control = rep(c(0,1), each = 50)), subset(unsat_gov_means, select=-c(ps_corp, state_control)), by="merge")
m2_unsat_gov_pred_df <- unique(m2_unsat_gov_pred_df)
#for m2_unsat_gov: ps_lib
m2b_unsat_gov_pred_df <- left_join(data.frame(merge = 1, outcome = "dissatisfied with government", type ="domestic", ps_lib = rep(seq(0, 1, length.out = 50), 2), state_control = rep(c(0,1), each = 50)), subset(unsat_gov_means, select=-c(ps_lib, state_control)), by="merge")
m2b_unsat_gov_pred_df <- unique(m2b_unsat_gov_pred_df)
#for m3_unsat_gov: ps_corp ps_corp_all_diff_neg
m3_unsat_gov_pred_df <- left_join(data.frame(merge = 1, outcome = "dissatisfied with government", type ="domestic", ps_corp = rep(seq(0, 1, length.out = 50), 3), ps_corp_all_max = rep(c(0,0.5,1), each = 50)), subset(unsat_gov_means, select=-c(ps_corp, ps_corp_all_diff_neg)), by="merge")
m3_unsat_gov_pred_df$ps_corp_all_diff_neg <- pmax(0, m3_unsat_gov_pred_df$ps_corp_all_max - m3_unsat_gov_pred_df$ps_corp)
m3_unsat_gov_pred_df <- unique(m3_unsat_gov_pred_df)
#for m4_unsat_gov ps_corp tek_ps_corp_diff_neg
m4_unsat_gov_pred_df <- left_join(data.frame(merge = 1, outcome = "dissatisfied with government", type ="TEK peers", ps_corp = rep(seq(0, 1, length.out = 50), 3), tek_ps_corp_all_max = rep(c(0,0.5,1), each = 50)), subset(unsat_gov_means, select=-c(ps_corp, ps_corp_all_diff_neg)), by="merge")
m4_unsat_gov_pred_df$tek_ps_corp_diff_neg <- pmax(0, m4_unsat_gov_pred_df$tek_ps_corp_all_max - m4_unsat_gov_pred_df$ps_corp)
m4_unsat_gov_pred_df <- unique(m4_unsat_gov_pred_df)

######2.3 predictions######
#for m1_unsat_gov: included_nc
m1_unsat_gov_pred_fun <- function(dat, inds)
{
  dat =  subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50)[inds,]
  m1_unsat_gov_temp <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = dat)
  pred.b <- predict(m1_unsat_gov_temp, newdata = m1_unsat_gov_pred_df, type ="response", re.form = NA)
}
m1_unsat_gov_pred <- boot(data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50), statistic = m1_unsat_gov_pred_fun, 100)
m1_unsat_gov_pred_tidy <- tidy(m1_unsat_gov_pred)
m1_unsat_gov_pred_df2 <- cbind(m1_unsat_gov_pred_df, m1_unsat_gov_pred_tidy)
#for m2_unsat_gov: ps_corp
m2_unsat_gov_pred_fun <- function(dat, inds)
{
  dat =  subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50)[inds,]
  m2_unsat_gov_temp <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = dat)
  pred.b <- predict(m2_unsat_gov_temp, newdata = m2_unsat_gov_pred_df, type ="response", re.form = NA)
}
m2_unsat_gov_pred <- boot(data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50), statistic = m2_unsat_gov_pred_fun, 100)
m2_unsat_gov_pred_tidy <- tidy(m2_unsat_gov_pred)
m2_unsat_gov_pred_df2 <- cbind(m2_unsat_gov_pred_df, m2_unsat_gov_pred_tidy)
#for m2_unsat_gov: ps_lib
m2b_unsat_gov_pred_fun <- function(dat, inds)
{
  dat =  subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50)[inds,]
  m2b_unsat_gov_temp <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = dat)
  pred.b <- predict(m2b_unsat_gov_temp, newdata = m2b_unsat_gov_pred_df, type ="response", re.form = NA)
}
m2b_unsat_gov_pred <- boot(data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50), statistic = m2b_unsat_gov_pred_fun, 100)
m2b_unsat_gov_pred_tidy <- tidy(m2b_unsat_gov_pred)
m2b_unsat_gov_pred_df2 <- cbind(m2b_unsat_gov_pred_df, m2b_unsat_gov_pred_tidy)
#for m3_unsat_gov: ps_corp ps_corp_all_diff_neg
m3_unsat_gov_pred_fun <- function(dat, inds)
{
  dat =  subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50)[inds,]
  m3_unsat_gov_temp <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = dat)
  pred.b <- predict(m3_unsat_gov_temp, newdata = m3_unsat_gov_pred_df, type ="response", re.form = NA)
}
m3_unsat_gov_pred <- boot(data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50), statistic = m3_unsat_gov_pred_fun, 100)
m3_unsat_gov_pred_tidy <- tidy(m3_unsat_gov_pred)
m3_unsat_gov_pred_df2 <- cbind(m3_unsat_gov_pred_df, m3_unsat_gov_pred_tidy)
#for m4_unsat_gov ps_corp tek_ps_corp_diff_neg
m4_unsat_gov_pred_fun <- function(dat, inds)
{
  dat =  subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50)[inds,]
  m4_unsat_gov_temp <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = dat)
  pred.b <- predict(m4_unsat_gov_temp, newdata = m4_unsat_gov_pred_df, type ="response", re.form = NA)
}
m4_unsat_gov_pred <- boot(data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50), statistic = m4_unsat_gov_pred_fun, 100)
m4_unsat_gov_pred_tidy <- tidy(m4_unsat_gov_pred)
m4_unsat_gov_pred_df2 <- cbind(m4_unsat_gov_pred_df, m4_unsat_gov_pred_tidy)

######2.4 plots######
#for m1_unsat_gov and m2_unsat_gov: included_nc ps_corp
m1_unsat_gov_pred_df3 <- m1_unsat_gov_pred_df2
m1_unsat_gov_pred_df3$group_type <- ifelse(m1_unsat_gov_pred_df3$group_type == "excluded NC", "non-core\n(excluded)", m1_unsat_gov_pred_df3$group_type)
m1_unsat_gov_pred_df3$group_type <- ifelse(m1_unsat_gov_pred_df3$group_type == "included NC", "non-core\n(de-facto included)", m1_unsat_gov_pred_df3$group_type)
m2_unsat_gov_pred_df2_corp <- subset(m2_unsat_gov_pred_df2, state_control == 0 & ps_corp == 1)
m2_unsat_gov_pred_df2_corp$group_type <- "non-core\n(corporate PSI)"
m2b_unsat_gov_pred_df2_lib <- subset(m2b_unsat_gov_pred_df2, state_control == 0 & ps_lib == 1)
m2b_unsat_gov_pred_df2_lib$group_type <- "non-core\n(liberal PSI)"
m12_comparison_unsat_gov_pred_df2 <- rbind.fill(m1_unsat_gov_pred_df3, m2_unsat_gov_pred_df2_corp, m2b_unsat_gov_pred_df2_lib)
m12_comparison_unsat_gov_pred_df2 <- subset(m12_comparison_unsat_gov_pred_df2, group_type != "core")
m12_comparison_unsat_gov_pred_df2$group_type <- as.factor(m12_comparison_unsat_gov_pred_df2$group_type)
m12_comparison_unsat_gov_pred_df2$group_type = factor(m12_comparison_unsat_gov_pred_df2$group_type,levels(m12_comparison_unsat_gov_pred_df2$group_type)[c(3,2,1,4)])
m12_comparison_unsat_gov_pred_plot <- ggplot(data = m12_comparison_unsat_gov_pred_df2) + 
  geom_point(aes(x = group_type, y = statistic)) +
  geom_errorbar(aes(x = group_type, ymin = statistic - 1.645 * std.error, ymax = statistic + 1.645 * std.error),width=0.1) +
  theme_bw() + theme(text=element_text(family='Times New Roman'), legend.position="bottom", plot.title = element_text(size=10), axis.title = element_text(size=10), axis.text = element_text(size=10)) + scale_y_continuous(labels = scales::percent_format(accuracy = 0.01)) +
  xlab('group type') + ylab('predicted probability\nof government\ndissatisfaction')
#for m3_unsat_gov: ps_corp ps_corp_all_diff_neg
m3_unsat_gov_pred_plot <- ggplot(data = m3_unsat_gov_pred_df2) + geom_line(aes(y=statistic, x=ps_corp, group = factor(ps_corp_all_max), colour = factor(ps_corp_all_max), linetype = factor(ps_corp_all_max))) + 
  geom_ribbon(alpha = 0.25, aes(x=ps_corp, ymin = statistic - 1.645 * std.error, ymax = statistic + 1.645 * std.error, group = factor(ps_corp_all_max), fill = factor(ps_corp_all_max))) +
  theme_bw() + theme(text=element_text(family='Times New Roman'), legend.position="bottom", plot.title = element_text(size=10), axis.title = element_text(size=10), axis.text = element_text(size=10)) + scale_y_continuous(labels = scales::percent_format(accuracy = 0.01)) +
  xlab('corporate PSI (group)') + ylab('predicted probability of\ngovernment dissatisfaction') +
  scale_colour_manual(values = c("#fee8c8","#fdbb84","#e34a33"), name="avg. corporate PSI (reference group)") +
  scale_fill_manual(values = c("#fee8c8","#fdbb84","#e34a33"), name="avg. corporate PSI (reference group)")+
  scale_linetype_manual(values = c(1,2,3), name="avg. corporate PSI (reference group)")
#for m4_unsat_gov ps_corp tek_ps_corp_diff_neg
m4_unsat_gov_pred_plot <- ggplot(data = m4_unsat_gov_pred_df2) + geom_line(aes(y=statistic, x=ps_corp, group = factor(tek_ps_corp_all_max), colour = factor(tek_ps_corp_all_max), linetype = factor(tek_ps_corp_all_max))) + 
  geom_ribbon(alpha = 0.25, aes(x=ps_corp, ymin = statistic - 1.645 * std.error, ymax = statistic + 1.645 * std.error, group = factor(tek_ps_corp_all_max), fill = factor(tek_ps_corp_all_max))) +
  theme_bw() + theme(text=element_text(family='Times New Roman'), legend.position="bottom", plot.title = element_text(size=10), axis.title = element_text(size=10), axis.text = element_text(size=10)) + scale_y_continuous(labels = scales::percent_format(accuracy = 0.01)) +
  xlab('corporate PSI (group)') + ylab('predicted probability of\ngovernment dissatisfaction') +
  scale_colour_manual(values = c("#fee8c8","#fdbb84","#e34a33"), name="avg. corporate PSI (reference group)") +
  scale_fill_manual(values = c("#fee8c8","#fdbb84","#e34a33"), name="avg. corporate PSI (reference group)")+
  scale_linetype_manual(values = c(1,2,3), name="avg. corporate PSI (reference group)")


####################################################################################
#########STEP 3: Predicted probability plots (feeling discriminated)#########
####################################################################################

######3.1 mean/median/modal values, all/included groups, used for predicted probabilities######
disc_grp_means <- subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50)
disc_grp_means <- subset(disc_grp_means, state_control == 0)
disc_grp_means <- disc_grp_means[c("survey","cowcode","cowcode_year","gwgroupid_year","state_control","included_nc","ps_corp","ps_corp_all_diff_neg","tek_ps_corp_diff_neg","ps_lib",controls_gc,controls_i)]
for (i in (categorical_vars)) {
  step1 <- paste("disc_grp_means$",i,"<- median(disc_grp_means$",i,")",sep="")
  eval(parse(text=c(step1)))
}
for (i in (nominal_vars)) {
  step1 <- paste("disc_grp_means$",i,"<- mlv(disc_grp_means$",i,")",sep="")
  eval(parse(text=c(step1)))
}
for (i in (numeric_vars)) {
  step1 <- paste("disc_grp_means$",i,"<- median(disc_grp_means$",i,")",sep="")
  eval(parse(text=c(step1)))
}
disc_grp_means <- unique(disc_grp_means)
disc_grp_means$merge <- 1

######3.2 dataframe: scenarios######
#for m1_unsat_gov: included_nc
m1_disc_grp_pred_df <- left_join(data.frame(merge = 1, outcome = "dissatisfied with government", type ="domestic", included_nc = c(0, 0, 1), state_control = c(1,0,0), group_type = c("core","excluded NC","included NC")), subset(disc_grp_means, select=-c(included_nc, state_control)), by="merge")
m1_disc_grp_pred_df <- unique(m1_disc_grp_pred_df)
#for m3_disc_grp: ps_corp
m2_disc_grp_pred_df <- left_join(data.frame(merge = 1, outcome = "feeling discriminated", type ="domestic", ps_corp = rep(seq(0, 1, length.out = 50), 2), state_control = rep(c(0,1), each = 50)), subset(disc_grp_means, select=-c(ps_corp, state_control)), by="merge")
m2_disc_grp_pred_df <- unique(m2_disc_grp_pred_df)
#for m3_disc_grp: ps_lib
m2b_disc_grp_pred_df <- left_join(data.frame(merge = 1, outcome = "feeling discriminated", type ="domestic", ps_lib = rep(seq(0, 1, length.out = 50), 2), state_control = rep(c(0,1), each = 50)), subset(disc_grp_means, select=-c(ps_lib, state_control)), by="merge")
m2b_disc_grp_pred_df <- unique(m2b_disc_grp_pred_df)
#for m3_disc_grp: ps_corp ps_corp_all_diff_neg
m3_disc_grp_pred_df <- left_join(data.frame(merge = 1, outcome = "feeling discriminated", type ="domestic", ps_corp = rep(seq(0, 1, length.out = 50), 3), ps_corp_all_max = rep(c(0,0.5,1), each = 50)), subset(disc_grp_means, select=-c(ps_corp, ps_corp_all_diff_neg)), by="merge")
m3_disc_grp_pred_df$ps_corp_all_diff_neg <- pmax(0, m3_disc_grp_pred_df$ps_corp_all_max - m3_disc_grp_pred_df$ps_corp)
m3_disc_grp_pred_df <- unique(m3_disc_grp_pred_df)
#for m4_disc_grp ps_corp tek_ps_corp_diff_neg
m4_disc_grp_pred_df <- left_join(data.frame(merge = 1, outcome = "feeling discriminated", type ="TEK peers", ps_corp = rep(seq(0, 1, length.out = 50), 3), tek_ps_corp_all_max = rep(c(0,0.5,1), each = 50)), subset(disc_grp_means, select=-c(ps_corp, ps_corp_all_diff_neg)), by="merge")
m4_disc_grp_pred_df$tek_ps_corp_diff_neg <- pmax(0, m4_disc_grp_pred_df$tek_ps_corp_all_max - m4_disc_grp_pred_df$ps_corp)
m4_disc_grp_pred_df <- unique(m4_disc_grp_pred_df)

######3.3 predictions######
#for m1_unsat_gov: included_nc
m1_disc_grp_pred_fun <- function(dat, inds)
{
  dat =  subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50)[inds,]
  m1_disc_grp_temp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = dat)
  pred.b <- predict(m1_disc_grp_temp, newdata = m1_disc_grp_pred_df, type ="response", re.form = NA)
}
m1_disc_grp_pred <- boot(data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50), statistic = m1_disc_grp_pred_fun, 100)
m1_disc_grp_pred_tidy <- tidy(m1_disc_grp_pred)
m1_disc_grp_pred_df2 <- cbind(m1_disc_grp_pred_df, m1_disc_grp_pred_tidy)
#for m2_disc_grp: ps_corp
m2_disc_grp_pred_fun <- function(dat, inds)
{
  dat =  subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50)[inds,]
  m2_disc_grp_temp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = dat)
  pred.b <- predict(m2_disc_grp_temp, newdata = m2_disc_grp_pred_df, type ="response", re.form = NA)
}
m2_disc_grp_pred <- boot(data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50), statistic = m2_disc_grp_pred_fun, 100)
m2_disc_grp_pred_tidy <- tidy(m2_disc_grp_pred)
m2_disc_grp_pred_df2 <- cbind(m2_disc_grp_pred_df, m2_disc_grp_pred_tidy)
#for m2_disc_grp: ps_lib
m2b_disc_grp_pred_fun <- function(dat, inds)
{
  dat =  subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50)[inds,]
  m2b_disc_grp_temp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = dat)
  pred.b <- predict(m2b_disc_grp_temp, newdata = m2b_disc_grp_pred_df, type ="response", re.form = NA)
}
m2b_disc_grp_pred <- boot(data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50), statistic = m2b_disc_grp_pred_fun, 100)
m2b_disc_grp_pred_tidy <- tidy(m2b_disc_grp_pred)
m2b_disc_grp_pred_df2 <- cbind(m2b_disc_grp_pred_df, m2b_disc_grp_pred_tidy)
#for m3_disc_grp: ps_corp ps_corp_all_diff_neg
m3_disc_grp_pred_fun <- function(dat, inds)
{
  dat =  subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50)[inds,]
  m3_disc_grp_temp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = dat)
  pred.b <- predict(m3_disc_grp_temp, newdata = m3_disc_grp_pred_df, type ="response", re.form = NA)
}
m3_disc_grp_pred <- boot(data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50), statistic = m3_disc_grp_pred_fun, 100)
m3_disc_grp_pred_tidy <- tidy(m3_disc_grp_pred)
m3_disc_grp_pred_df2 <- cbind(m3_disc_grp_pred_df, m3_disc_grp_pred_tidy)
#for m4_disc_grp ps_corp tek_ps_corp_diff_neg
m4_disc_grp_pred_fun <- function(dat, inds)
{
  dat =  subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50)[inds,]
  m4_disc_grp_temp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = dat)
  pred.b <- predict(m4_disc_grp_temp, newdata = m4_disc_grp_pred_df, type ="response", re.form = NA)
}
m4_disc_grp_pred <- boot(data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50), statistic = m4_disc_grp_pred_fun, 100)
m4_disc_grp_pred_tidy <- tidy(m4_disc_grp_pred)
m4_disc_grp_pred_df2 <- cbind(m4_disc_grp_pred_df, m4_disc_grp_pred_tidy)

######3.4 plots######
#for m1_unsat_gov and m2_unsat_gov: included_nc ps_corp
m1_disc_grp_pred_df3 <- m1_disc_grp_pred_df2
m1_disc_grp_pred_df3$group_type <- ifelse(m1_disc_grp_pred_df3$group_type == "excluded NC", "non-core\n(excluded)", m1_disc_grp_pred_df3$group_type)
m1_disc_grp_pred_df3$group_type <- ifelse(m1_disc_grp_pred_df3$group_type == "included NC", "non-core\n(de-facto included)", m1_disc_grp_pred_df3$group_type)
m2_disc_grp_pred_df2_corp <- subset(m2_disc_grp_pred_df2, state_control == 0 & ps_corp == 1)
m2_disc_grp_pred_df2_corp$group_type <- "non-core\n(corporate PSI)"
m2b_disc_grp_pred_df2_lib <- subset(m2b_disc_grp_pred_df2, state_control == 0 & ps_lib == 1)
m2b_disc_grp_pred_df2_lib$group_type <- "non-core\n(liberal PSI)"
m12_comparison_disc_grp_pred_df2 <- rbind.fill(m1_disc_grp_pred_df3, m2_disc_grp_pred_df2_corp, m2b_disc_grp_pred_df2_lib)
m12_comparison_disc_grp_pred_df2 <- subset(m12_comparison_disc_grp_pred_df2, group_type != "core")
m12_comparison_disc_grp_pred_df2$group_type <- as.factor(m12_comparison_disc_grp_pred_df2$group_type)
m12_comparison_disc_grp_pred_df2$group_type = factor(m12_comparison_disc_grp_pred_df2$group_type,levels(m12_comparison_disc_grp_pred_df2$group_type)[c(3,2,1,4)])
m12_comparison_disc_grp_pred_plot <- ggplot(data = m12_comparison_disc_grp_pred_df2) + 
  geom_point(aes(x = group_type, y = statistic)) +
  geom_errorbar(aes(x = group_type, ymin = statistic - 1.645 * std.error, ymax = statistic + 1.645 * std.error),width=0.1) +
  theme_bw() + theme(text=element_text(family='Times New Roman'), legend.position="bottom", plot.title = element_text(size=10), axis.title = element_text(size=10), axis.text = element_text(size=10)) + scale_y_continuous(labels = scales::percent_format(accuracy = 0.01)) +
  xlab('group type') + ylab('predicted probability\nof feeling\ndiscriminated')
#for m3_disc_grp: ps_corp ps_corp_all_diff_neg
m3_disc_grp_pred_plot <- ggplot(data = m3_disc_grp_pred_df2) + geom_line(aes(y=statistic, x=ps_corp, group = factor(ps_corp_all_max), colour = factor(ps_corp_all_max), linetype = factor(ps_corp_all_max))) + 
  geom_ribbon(alpha = 0.33, aes(x=ps_corp, ymin = statistic - 1.645 * std.error, ymax = statistic + 1.645 * std.error, group = factor(ps_corp_all_max), fill = factor(ps_corp_all_max))) +
  theme_bw() + theme(text=element_text(family='Times New Roman'), legend.position="bottom", plot.title = element_text(size=10), axis.title = element_text(size=10), axis.text = element_text(size=10)) + scale_y_continuous(labels = scales::percent_format(accuracy = 0.01)) +
  xlab('corporate PSI (group)') + ylab('predicted probability of\nfeeling discriminated') +
  scale_colour_manual(values = c("#fee8c8","#fdbb84","#e34a33"), name="avg. corporate PSI (reference group)") +
  scale_fill_manual(values = c("#fee8c8","#fdbb84","#e34a33"), name="avg. corporate PSI (reference group)")+
  scale_linetype_manual(values = c(1,2,3), name="avg. corporate PSI (reference group)")
#for m4_disc_grp ps_corp tek_ps_corp_diff_neg
m4_disc_grp_pred_plot <- ggplot(data = m4_disc_grp_pred_df2) + geom_line(aes(y=statistic, x=ps_corp, group = factor(tek_ps_corp_all_max), colour = factor(tek_ps_corp_all_max), linetype = factor(tek_ps_corp_all_max))) + 
  geom_ribbon(alpha = 0.33, aes(x=ps_corp, ymin = statistic - 1.645 * std.error, ymax = statistic + 1.645 * std.error, group = factor(tek_ps_corp_all_max), fill = factor(tek_ps_corp_all_max))) +
  theme_bw() + theme(text=element_text(family='Times New Roman'), legend.position="bottom", plot.title = element_text(size=10), axis.title = element_text(size=10), axis.text = element_text(size=10)) + scale_y_continuous(labels = scales::percent_format(accuracy = 0.01)) +
  xlab('corporate PSI (group)') + ylab('predicted probability of\nfeeling discriminated') +
  scale_colour_manual(values = c("#fee8c8","#fdbb84","#e34a33"), name="avg. corporate PSI (reference group)") +
  scale_fill_manual(values = c("#fee8c8","#fdbb84","#e34a33"), name="avg. corporate PSI (reference group)")+
  scale_linetype_manual(values = c(1,2,3), name="avg. corporate PSI (reference group)")


####################################################################################
#########STEP 4: PLOTS#########
####################################################################################

###a) figure 3
figure3 <- grid.arrange(m12_comparison_unsat_gov_pred_plot + ggtitle("a) government dissatisfaction"), m12_comparison_disc_grp_pred_plot + ggtitle("b) feeling discriminated"), nrow=2, ncol=1)

###b) figure 4
figure4 <- grid_arrange_shared_legend(m3_unsat_gov_pred_plot + ggtitle("a) government dissatisfaction,\ndomestic reference group"), m3_disc_grp_pred_plot + ggtitle("b) feeling discriminated,\ndomestic reference group"), m4_unsat_gov_pred_plot + ggtitle("c) government dissatisfaction,\ninternational (TEK) reference group"), m4_disc_grp_pred_plot + ggtitle("d) feeling discriminated,\ninternational (TEK) reference group"), nrow=2, ncol=2)


